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We consider the simplest example of lattice frustration in the | -filled band, a one-dimensional chain with next-nearest 
neighbor interactions. For this zigzag ladder with electron-electron as well as electron-phonon interactions we present 
numerical results for ground state as well as thermodynamic properties. In this system the ground state bond distortion 
pattern is independent of electron-electron interaction strength. The spin gap from the ground state of the zigzag ladder 
increases with the degree of frustration. Unlike in one-dimension, where the spin-gap and charge ordering transitions can 
be distinct, we show that in the ladder they occur simultaneously. We discuss spin gap and charge ordering transitions 
in I -filled materials with one, two, or three dimensional crystal structures. We show empirically that regardless of 
dimensionality the occurrence of simultaneous or distinct charge and magnetic transitions can be correlated with the 
ground state bond distortion pattern. 
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1. Introduction 

The consequences of geometric lattice frustration on 
ground and excited state behavior of spin systems have been 
widely studied over the past two decades.''^ Even as many 
of the issues are still being debated for the complex two- and 
three-dimensional (2D and 3D) lattices, consensus on some 
of the simpler models have been reached. One such model 
frustrated system is the spin-| zigzag ladder, with nonzero 
antiferromagnetic exchanges Ji and J2 between nearest and 
next-nearest neighbor spins. The ground state here is a va- 
lence bond solid (VBS) for J2 - 0.57i,^'^ and is dimerized for 
J2 > 0.241 17], ^'^ where the excitations are spin solitons.^'^ 
Some theoretical results exist also for J2 > J\ ' 

Spin Hamiltonians are restricted to a charge per site p - \ 
(^-filled band) and are obtained in the limit of very large on- 
site Hubbard repulsion U between electrons or holes where 
there are no charge degrees of freedom. The literature on 
frustrated non-^-filled bands, where the charge degrees of 
freedom are non-vanishing even for f/ — > 00, remains rel- 
atively sparse. Much of the literature on non-^-filled band 
frustrated systems is for p - 5, or the ^-filled band, which 
is of interest both because under certain conditions it can be 
described within an effective i-filled band picture,'^ and 
because there exist many i-filled band frustrated materials, 
both organic and inorganic, that exhibit novel behavior includ- 
ing charge, spin and orbital-ordering, and superconductivity 
(see below). We have recently shown that there is a strong 
tendency to form local spin-singlets in the | -filled band with 
strong electron-electron interactions, in both one dimension 
(ID) and 2D, and especially in the presence of lattice frus- 
tration, which enhances quantum effects.'^- The frustration- 
driven singlet formation does not occur if only the frustrat- 
ing Coulomb interactions are included;'^ the inclusion of the 
frustrating electron-hopping integral, with or without the cor- 
responding Coulomb interactions is essential. Singlet forma- 
tion in the | -filled band is accompanied by charge-ordering 




Fig. 1. Zigzag ladder lattice and ground state configurations for density 
p = 5. (a) The ground state for (2/?! > 1.707 has uniform charges and bonds, 
(b) For ?2/f| < 1-707, the ground state is a PEC with coexisting charge and 
bond order. Filled (open) circles correspond to sites with large (small) charge 
density. The heavy line indicates the strongest (spin-singlet) bond. Double, 
single, and dotted lines indicate strong, intermediate, and weak bonds, re- 
spectively. 



(CO), leading to what we have termed as a paired-electron 
crystal (PEC), which consists of pairs of singly-occupied sites 
separated by pairs of vacant sites. The PEC is the ^-filled 
band equivalent of the VBS. The difference from the stan- 
dard VBS lies in the possibility that the PEC gives way to 
a paired-electron liquid under appropriate conditions, which 
might then lead to superconductivity. Clearly it is desirable to 
extend of our cuiTent work on the ground state of the PEC to 
excited states. 

In the present paper, we discuss theoretical results for the 
PEC that occurs in the simplest ^-filled frustrated lattice with 
interacting electrons. The lattice geometry here is the same 
as in the spin zigzag ladder (see Fig. 1(a)), with nearest and 
next nearest neighbor electron hopping matrix elements t\ and 
f2, respectively. In our previous work on the ground state of 
the ^-filled band zigzag ladder,'^ we had shown that in the 
presence of electron-lattice coupling (either nearest neighbor 
along the zigzag direction, or second neighbor, or both), and 
for f2/f 1 less than a critical value (f2/fi)c (see below) there oc- 
curs a coupled bond-charge density wave (BCDW), as shown 
in Fig. 1(b). Alternate rungs of the zigzag ladder are occu- 
pied by singlet spin-coupled charge-rich sites, and charge- 
poor sites occupy the other set of rungs. The CO pattern here 
characterizes the BCDW as a PEC. We discuss below the na- 
ture of the ground state in this system for the complete param- 
eter ranee hlu < (hiti).. We also oresent numerical results 
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for spin excitations and thermodynamics for the | -filled band 
zigzag ladder within the interacting electron Hamiltonian and 
discuss applications of the model to real systems. 

Although our actual work is limited to the zigzag ladder, 
our analysis leads to an empirical understanding of a perplex- 
ing observation in the ^-filled band systems with interact- 
ing electrons in general. In many such systems, irrespective 
of dimensionality, there often are two transitions, (i) a high 
temperature metal-insulator (MI) transition at Tmi that is ac- 
companied by CO without any perceptible effect on the mag- 
netic behavior, and (ii) a second insulator-insulator transition 
at a lower temperature Tsg where a spin gap (SG), seen in 
the magnetic susceptibility, occurs. In a second class of sys- 
tems with nearly identical chemical constituents there occurs 
however a single transition with Tmi = Tsg- We point out a 
well-defined correlation between this "two-versus-one" tran- 
sition and the pattern of the bond distortion in the spin-gapped 
phase. Furthermore, many interacting i-fiUed band materi- 
als are superconducting, often under pressure. There appears 
to exist a correlation also between superconductivity and the 
bond distortion pattern at the lowest temperatures. 

The organization of the paper is as follows. In Section 2 we 
introduce our Hamiltonian for the ^-filled band zigzag ladder. 
In Section 3 we present a brief summary of the earlier results 
for the ID limit of this model, ^^"^^ and for the ground state 
of the zigzag ladder."^ Two different bond distortion patterns, 
accompanied however with the same charge distortion pattern 
are possible in the ID limit. We point out the two bond dis- 
tortion patterns correspond to two different mappings of the 
^-filled band to effective i-fiUed bands. In the ID limit both 
mappings are valid. In contrast, only one mapping is appli- 
cable to the zigzag ladder, and there is a single unambigu- 
ous PEC here. The results of our numerical calculations for 
the zigzag ladder are presented in Section 4, where we dis- 
cuss both ground state and temperature-dependent behavior. 
In Section 5 we discuss the implications of our results for real 
I -filled materials. Appendix A contains details of the numer- 
ical method, based on a matrix-product state (MPS) repre- 
sentation, used for our large-lattice calculations. Appendix B 
contains details of the finite-size scaling of our ground-state 
calculations. 

2. Theoretical model and parameter space 

The Hamiltonian for the zigzag ladder is 

i i 

- t2 ^(1 + a2Ai,i+2)Bi,i+2 + ^K2 ^li+2 
i i 

+ ^{Uni^riii + ViHirii+i + ¥2^1^+2)- (1) 

In Eq. 1 Bi j - E(7-(cj,^c,,o- + H.c.) is the kinetic energy op- 
erator for the bond between sites i and /, where c! creates 

an electron of spin cr on site i. riia- = c]^Ci^a- is the density 
operator, and n,- = n,-f + ^4. ti and t2 are hopping integrals 
along the zigzag rung direction and the rail direction, respec- 
tively, as shown in Fig. 1. The lattice may also be viewed as 
a ID chain with nearest neighbor hopping f 1 and frustrating 



below, however, we will imply t2 - and V2 - 0. We give 
energies in units of ti and fix fi = l. A,_y is the deformation of 
the bond between sites i and j; a\ and 02 are the inter-site 
electron-phonon (e-p) coupling constants with corresponding 
spring constants Ki and K2, which for simplicity we choose 
to be identical in the t\ and t2 directions (ori = a2 = a and 
Ki = K2 = K). For all results below we fix Ki = K2 - 2. 
We have omitted intra-site e-p coupling terms in Eq. 1, as 
apart from changing the strength of the charge order, they 
do not greatly change the thermodynamic properties. U is 
the onsite Coulomb interaction and Vi and V2 are the nearest- 
neighbor Coulomb interactions for ti and t2 bonds. We will 
consider only the case V] - V2 - V. In the ID limit, the 
ground state is a Wigner crystal (WC) for Vi > V^, where 
= 2\ti\ for f/ -> 00 and is larger for finite U. Our discus- 
sions of the ID limit are for Vi < V^. All calculations use 
periodic boundary conditions. 

3. Summary of earlier results. 

Given the complexity of our results, it is useful to have a 
brief summary of our earlier results. While the numerical re- 
sults for the ID limit and the ground state of the zigzag ladder 
have been presented before, the mappings of the ground states 
to the effective ^-filled band model that we discuss below give 
new insight. 

3.1 ID limit 

For moderate to strong electron-electron (e-e) interactions 
but V\ < Vp the ^-filled band is bond- or charge-dimerized 
at T<Tmi. In either case the ground state enters a spin-Peierls 
(SP) phase^^'^-' for T<Tsg. We show in Fig. 2(a) a schematic 
of the bond-dimerized phase, with site charge densities 0.5, 
strong intradimer bonds and weak interdimer bonds. The 
dimer unit cells (the boxes in the figure) containing a single 
electron can be thought of as effective sites, which leads to the 
effective | -filled band description in this case.''* The descrip- 
tion of the SP state is then as shown in Fig. 2(b), with alter- 
nating /nfer-dimer bonds. The overall pattern of bond strength 
here is ■ • • Strong- Weak-Strong-Weak' • • • (SWSW). 
UnUke the true ^ -filled band, however, the charge degrees 
of freedom internal to the dimer unit cell are relevant in 
the i -filled band, and the charge distribution is as shown in 
the Fig. 2(b). The coexisting charge order pattern is written 
as •••llOO---, where '1' ('0') denotes site charge density 
of 0.5-1-5 (0.5-6). The CO ampHtude 6 varies with e-e and 
e-p interaction strengths. Nearest-neighbor spin-singlet cou- 
pling between the electrons on charge-rich '1' sites that are 
linked by the W' inter-dimer bond gives the spin gap of this 
SP phase. This coexisting broken symmetry state has been 
termed a BCDW,'^'^''^'* and is the simplest example of the 
PEC found more generally in | -filled systems beyond the ID 
limit. Finite temperature phase transition to the SWSW' 
phase, starting from the uniform phase and through the dimer 
phase of Fig. 2(a), has been demonstrated within an extended 
Hubbard model that incorporated interchain interactions at a 
mean-field level.^^ 

The pattern of bond order modulation in the ID PEC is 
not unique but depends on the strength of e-e interactions.'^ 
The bond-order pattern that dominates for relatively weaker 
e-e interactions, or for relatively strong e-p interaction has the 



J. Phys. Soc. Jpn. 



FULL PAPERS 



(a) [y^^sTsV^sTsKs^ (c) t=|— 0— 1=|— o— 

(b) \W\- \W] rQ^=flo] (d)2— — 2— — 

Fig. 2. (a) Bond-dimerized phase with uniform charge density. Boxes de- 
note dimer units, (b) SP state evolving from (a), with bond pattern SWSW. 
Strong 'S' dimer bonds are represented by boxes enclosing two sites. Zeroes 
indicate sites with charge density 0.5-5 while spins indicate sites with charge 
density 0.5+5. Double dotted bonds are stronger than single-dotted bonds, 
but weaJcer than the dimer bonds, (c) Bond pattern SMWM found in the case 
of relatively weaJcer e-e interactions. Here double lines indicate the strongest 
'S' bonds, single lines medium strength 'M' bonds, and dotted lines weak 
'W bonds, (d) Effective ^ -filled band model for (c). 



portantly, in this bond pattern the coexisting charge modula- 
tion pattern is again ■ ■ ■ 1 100- ■ ■ but the spin singlets now co- 
incide with the strongest 'S' mfra-dimer bonds (Fig. 2(c)). As 
shown in Fig. 2(d), the corresponding effective i-filled band is 
now different; pairs of sites with large (small) charge densities 
are now mapped onto effective sites with charge density 2 (0). 
The bonds between the effective sites are now equivalent, and 
the effective state is now a ^-filled band .s/fe-diagonal CDW, 
which is known to have a single transition where charge and 
magnetic gaps open simultaneously. Further, as the strongest 
bond coincides with the location of the spin-singlet, this state 
is also expected to have a larger SG and transition temperature 
than TsG in the SWSW case. A simple physical picture for 
two-versus-one transition is thus obtained from the ID study: 
if the strongest bonds are between a charge-rich and a charge- 
poor site, the spin-singlet formed are inter-dimer and there 
will occur two transitions; if, however, the strongest bond is 
a 1-1 bond, the spin singlet is located on an intra-dimer bond 
and there is a single transition. 

3.2 Zigzag ladder 

The energy dispersion relation in the zigzag ladder in the 
noninteracting limit is a sum of two terms, 

E{k) = -2f 1 cos(^) - 2f2 cos(2^). (2) 

Here q refers to a wavenumber along the fj direction, view- 
ing the ladder as a ID chain with second-neighbor interac- 
tions. The topology of the bandstructure changes at f2/fi - 
(t2/h)c = (2+ V2)/2 = 1.707 ■■■; for f2/fi < (f2/fi)c- the 
Fermi surface consists of two points aX q — kf - j, while for 
hiti > (t2/ti)c there are four such points. In the presence of 
e-p coupling the ground state is then unstable'** to a Peierls 
distorted state for fa/fi < (f2/fi)c. As in ID, the ground state 
in the distorted region again has ■ ■ ■ 1 100 ■ ■ ■ CO and a spin 

gap- 
Importantly, for the p = 5 zigzag ladder, there is no WC 
state with a charge order distinct from the ■ ■ ■ 1100 ■ ■ ■ charge 
order found in the PEC (for V] ~ ¥2)- In the zigzag ladder, 
the WC charge pattern ■ ■ ■ 1010 ■ ■ ■ CO can be placed along 
the two f2 directions in two different ways, both of which lead 
to the same PEC state shown in Fig. 1(b). Placing the WC CO 
pattern ■ ■ ■ 1010 ■ ■ ■ along the zigzag direction would place all 
charge on a single t2 chain and is stable only in the limit Vi » 
V2- Thus only the PEC of Fig. 1(b) is stable for realistic Vi ~ 
V2. As we show below, this has important consequences for 
both the number of transitions expected in the zigzag ladder 
and the nature of the spinon excitations. 




Fig. 3. (color online) (a) E,, the per site spring constant energy term in 
Eq. 1, as a function of b/fi for a 256 site ladder with U = V = 0. Electron- 
phonon parameters are ai = 1.6, Qr2 = (circles); ai = 0, Q'2 = 1.6 (squares); 
and a I = 1.6, ai = \.(> (diamonds), (b) Cooperative enhancement Ai?, as a 
function of tilti (see text). The inset shows the small (2/?! region in more 
detail. 



4. Numerical Results for the zigzag ladder 

We use two different numerical methods to solve Eq. 1 . For 
the ground state order parameters and the spin gaps in the in- 
teracting case we use a new variational quantum Monte Carlo 
(QMC) method using a MPS basis. For quasi-lD systems, 
this MPS-QMC method provides accuracy similar to that 
of the Density Matrix RenormaUzation Group (DMRG)^^'^^ 
method. Similar methods using MPS representations have 
predominantly been used to study quantum spin systems. Our 
results here show however that they may be successfully ap- 
plied to electronic models as well. Details on the application 
of this method to Hubbard-type models are discussed in Ap- 
pendix A. For finite temperatures our calculations are for zero 
e-p interactions {a\ - 02 - 0). This is because information 
about the tendency to distortion exists in the wavefunction 
even without inclusion of explicit e-p interactions.^'' We use 
here the standard determinantal QMC method.^" While the 
lowest temperatures achievable by this method are limited by 
the Fermion sign problem, in the present system for density 
P - \, inverse temperatures ofj8 * 10 - 16 are reachable for 
parameters U ~ 4 and V = 0. 

4.1 Charge order, bond periodicity, and spin gap 

We first consider the solution of Eq. 1 in the limit of zero 
e-e interactions ([/ = V = 0). Because Eq. 1 includes e-p 
coupling in both the fi and f2 directions, the variation of the 
lattice distortion amplitude with the ratio f2/fi is nontrivial. 
The total spring constant energy Eg, the sum of Ki and K2 
terms in Eq. 1, is a convenient measure of the strength of the 
lattice distortion. Fig. 3(a) shows as a function of f2/fi for 
a 256 site ladder Three different choices for e-p couplings are 
shown: (i) ai - a, 02 - 0; (ii) = 0, a'2 = a; and (iii) 
ai = a2 - a. In the first case with e-p coupling only along 
the zigzag direction, the lattice distortion is strongest when 
t2 = 0, and vanishes continuously at {t2/t\)c- In the second 
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Fig. 4. Order parameters versus f2/?i for the zigzag ladder with U = 6, 

V = 1, and a = 1.6. Results are finite-size scaled from 20, 28, 36, and 60 site 
lattices, (a) 4kf component of the lattice distortion, and (b) spring energy 
per site, Es, (c) charge disproportionation An, and (d) singlet-triplet gap Ast. 
Open symbols are for V = Vi = V2, and filled symbols at t2/ti = are for 

V = Vi, ¥2 = (see text). 



case with e-p coupling only along the ?2 bonds, the strength 
of the distortion increases with f2- Putting these two effects to- 
gether, the strength of the lattice distortion in the general case 
with ori and a2 nonzero shows a minimum at an intermediate 
value of t2/ti, as seen in Fig. 3. The precise tj/ti value of this 
minimum will depend on the specific choices for ai and ai. 

Fig. 3(a) also shows that the ?2 and t\ lattice distortions act 
cooperatively-the total lattice distortion for both ai and ck2 
nonzero is considerably stronger than the sum of the inde- 
pendent distortions. In Fig. 3(b), we plot a quantitative mea- 
sure of the cooperative enhancement, AEg = Es(ai - a2 - 
a) - {Es{a\ = a,a2 = 0) + Esiai = 0, ai = a)), as a func- 
tion of t2/h. AEs increases continuously with t2/h, and as 
shown in the inset is nonzero even for small values of 12/ ti. 
As we will discuss further in Section 5, this cooperative eff'ect 
also has potentially important consequences for materials — if 
both ai and a2 are nonzero, near {t2lh)c the lattice distortion 
changes abruptly. 

A key difference between the p = ^ zigzag ladder and the 
single chain is that in the ladder the bond pattern of the PEC 
remains SMWM regardless of the strength of e-e interactions. 
In ID the displacement of site j from equiUbrium, Uj {Ajj+\ = 
Uj+\ - Uj in Eq. 1), can be written as'' 

Uj = uo[r2 cos(2fcF j - 62) + rn cos{Ak^ j - 64)], (3) 

where r2 and are relative components of the period-4 2kF 
and period-2 4kF lattice distortions and mq is an overall ampli- 
tude. ^2 and ^4 are normaUzed such that r2 + r^^ = 1. The 
SMWM bond pattern corresponds to ra, < 0.41 while the 
SWSW pattern''^ has > 0.41. In both cases 6*2 = | and 
04 = 0. Note that Ajj+2 in Eq. 1 is an independent lattice dis- 
tortion of the t2 bonds and has no counterpart in the ID model. 

While the PEC occurs in the thermodynamic limit for 
t2lt\ < {t2lh)c and a = 0"^, in finite lattices a finite e-p 
coupling is required to observe the broken-symmetry ground 
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Fig. 5. Dependence of order parameters on V. Results are finite-size scaled 
from 20, 28, 36, and 60 site lattices. Quantities plotted in panels (a)-(d) are 
as defined in Fig. 4. Here t2/ti = 1.5, U = 6, and a = 1.6. 



larger than the minimum required for the ground and triplet 
states to be Peierls distorted. The strength of the lattice dis- 
tortion depends on the values of a and yS and therefore results 
here should not be directly compared to experimental values. 

In Fig. 4 we plot the finite-size scaled values of several or- 
der parameters determined self-consistently versus t2/ti for 
a = 1.6, U = 6, and V = I. Exact and MPS-QMC calcu- 
lations were performed for 20, 28, 36, and 60 site lattices. 
Details of the finite-size extrapolation are given in Appendix 
B. 

Fig. 4(a) shows the 4k/r lattice distortion strength and 
Fig. 4(b) the spring constant energy per site (including both ti 
and ?2 lattice distortions). Fig. 4(c) shows the charge dispro- 
portionation An, defined as the difference between the charge 
density on charge-rich and charge-poor sites. The singlet- 
triplet gap defined as Ast = E{S = 1) - E{S = 0) is plot- 
ted in Fig. 4(d). Note that two different 'ID' limits are shown 
in Fig. 4 at /2/?i = 0: either including the second neighbor 
Coulomb interaction {V2 = V, open symbols), or with only 
nearest-neighbor Coulomb interactions (¥2 = 0, filled sym- 
bols). 

Focusing first on the ^2/^1 = Umit, the bond pattern is 
SMWM in the case where V2 is nonzero, as is found in the 
ID weakly correlated band.'' In the more traditional ID with 
no second-neighbor V2, the bond distortion pattern is SWSW' 
for U - 6, V = 1 with ^4 > 0.41 (filled circle). For equal 
e-p coupling, the lattice distortion energy. An, and Ast are all 
stronger when the bond pattern is SMWM. Away from the ID 
limit, r4 increases with increasing 12/ ti but is always less than 
the 0.41 that would be necessary to reach the SWSW' bond 
pattern. The amplitude of the lattice distortion measured by 
the spring constant energy in Fig. 4(b) shows a minima at in- 
termediate t2/ti as in the non-interacting case. Surprisingly, 
An decreases as the lattice distortion becomes stronger-this 
is one significant difference from the ID chain BCDW state 
where An follows the strength of the bond distortion. Refer- 
ence'^ showed that the spin gap in the PEC state of the zigzag 
ladder is larger than the gap in the single chain having the 

camp An Ac Pie 4^'r^^ chnwc thp ear* in thp laHHpr inrrpacpc 
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with t2lt\, and is largest near (f2/fi)c- 

The U interaction weakens the PEC in the zigzag ladder 
(not shown). This decrease is, however, weaker than in the ID 
chain. In Fig. 5 we show the results of varying V while keep- 
ing other parameters fe/^i = 1.5, t/ = 6, and a = 1.6) fixed. 
Even at large V {V > f//2), as expected we did not find any 
transition to a distinct WC state; in all cases the CO pattern is 

• • • 1 100 • • • along the zigzag direction. Contrary to what oc- 
curs in the ID limit, the 4kf component of the bond distortion 
actually decreases with increasing V. While in ID V destabi- 
lizes the • • • 1 100 • • • CO, in the zigzag ladder, this effect in the 
ladder along the t\ direction is countered by V2, which prefers 

• • ■ 1010 ■ ■ ■ order along the t2 directions. A similar effect is 
found in the | -filled PEC state in the 2D anisotropic triangu- 
lar lattice — ^there also V strengthens the • • • 1 100 • • • CO pro- 
vided V is not too large. The most striking result of Fig. 5 
is the very large spin gap that is obtained when both 12! ti and 
V are moderately large. We will return to this point later in 
Section 5. 

4. 2 Thermodynamics and spinon binding 

To understand the thermodynamics of the zigzag ladder we 
present our results from complementary calculations of (i) 
wavenumber-dependent susceptibiUties, and (ii) the nature of 
higher spin states. Finite-temperature calculations of suscep- 
tibilities here are done within the static undistorted lattice. 

4.2.7 Temperature dependence of susceptibiUties 

We calculate wavenumber dependent charge OfpC^)), spin 
(^o-(^)), and bond susceptibiUties ixBiq)), defined as 

^-(^) = i Z r dT{0]{T)Oim. (4) 

In Eq. ^,Cf.= nj,t+nj,i, CF. = n^.t -n^u, and = Bjj+u for 
charge, spin, and bond order susceptibilities, respectively, p is 
the inverse temperature in units of ti. To facilitate comparison 
with ID, the q is again taken to be one dimensional as in Eq. 2. 
The presence and periodicity of charge and bond order can 
be determined from divergences of the charge or bond-order 
susceptibility as -> 00 Within determinantal QMC methods 
finite V leads to a significantly worse Fermion sign problem 
and hence we present results only for nonzero U but V = 0. 
However, as shown in Fig. 5(a), V does not change the pattern 
of the bond distortion. We therefore expect that our results 
here are representative for arbitrary U and V. 

Fig. 6 shows the bond, charge, and spin susceptibilities as 
a function of q and temperature for a = 36 site ladder with 
U = 4, y = 0, and two representative values of /2/?i, 0.4 
(Fig. 6(a)-(c)) and 1.4 (Fig. 6(d)-(f)). At low temperature, the 
bond, charge, and spin susceptibilities all peak at 2kf, consis- 
tent with period-four order of charges and bonds as expected. 
The 2kp peak in the spin susceptibiUty corresponds to the 
short-range AFM spin order found in the ID | -filled chain. As 
seen in Fig. 6(c) and (f), the spin susceptibility converges to a 
finite value as g — > at low temperatures. This is consistent 
with the expectation that no SG exists in the ladder in the limit 
of zero e-p coupling. Several differences can be noted when 
comparing results for small and large 12/ h . First, the Ikp bond 
and charge response as T— > is stronger for small t2/t]. This 
reflects the larger amnlinide distortion (An and bond order^ 



found in the ID weakly-correlated limit (see Fig. 4(b)-(c)). 
The 2kF spin response becomes weaker with increasing 12/ ti, 
and Xa-il) is reduced for small q. These changes reflect the 
increasing strength of the spin-singlet bond along the ti direc- 
tion with increasing t2/ti, which in turn leads to the increase 
in AsT- For larger t2/h, a broad plateau appears at ^ « 3n/4 
in the charge and spin susceptibilities. This plateau is however 
non-divergent — while it is significant at high and intermedi- 
ate T, at low temperatures the 2kF response becomes stronger. 
We will show below that this plateau reflects the binding of 
spinon excitations. 

Summarizing, the above numerical results show that in the 
zigzag ladder, no separate high-temperature ordering is ex- 
pected; instead the ladder is metallic at high temperature and 
as temperature decreases a single transition to the PEC state 
with SG takes place. This is in contrast to what happens in 
ID.22 

4.2.2 Spinon excitations 

The structure of excitations out of the ground state pro- 
vides an alternate way to understand the thermodynamics of a 
strongly-correlated system. In the ID limit, flipping one spin 
in the PEC results in two spinous which are unbound and 
hence separated by a distance of N/2 sites on the lattice. We 
have performed self-consistent calculations within Eq. 1 of 
excited spin triplet (5 = 1) excitations to detect and study 
such spinon excitations in the zigzag ladder. Our calculations 
are for moderately large e-p interactions so that the widths 
of the charge-spin solitons are relatively narrow. This is nec- 
essary in order to prevent overlaps between the solitons. In 
Fig. 7 we show the charge and spin densities in the 5 = 1 state 
of a 36 site zigzag ladder for small and large t2/t\. Spinon ex- 
citations are identified as defects in the ■ ■ ■ 1100 ■ ■ • charge 
order, and also from their large local spin densities. 

For small fi/fi (Fig. 7(a)), flipping a single spin results in 
two separated spinous as in the ID chain. The repulsive in- 
teraction between spinous results in their separating to oppo- 
site positions on the periodic lattice. As indicated in the fig- 
ure, each spinon occupies two lattice sites, and there occur 
both charge and spin modulations. The charge density at each 
spinon site is 0.5; the spin densities are also equal on the two 
sites. The PEC charge and bond distortions on either side of 
a spinon are out of phase with each other by two lattice sites. 
The charge densities on the two sites at the center of the are 
0.5, 0.5 (see Fig. 7(a)). With increasing t2/ti, the lattice sepa- 
ration between the spinous decreases, indicating binding. For 
the e-p coupling of Fig. 7, at approximately f2/f2 ~ 1.0 the 
bound spinous form a single excitation that occupies 4 lat- 
tice sites with approximately uniform charge density of 0.5 
on each site (Fig. 7(b)). For spin states higher than 1 (not 
shown here), we found that spinous are bound in pairs; for 
example in the 5 =2 state there are two of the defects shown 
in Fig. 7(b), separated by the maximum possible lattice spac- 
ing. Spinon binding is usually associated with an increase in 
the singlet-triplet gap, and can thus explain the observation 
the(Fig. 4(d)) that for fixed V, Ast increases with increasing 
f2/fi, even though An decreases at the same time. 

The thermodynamic behavior is understood only when the 
complementary calculations of Fig. 6 and 7 are taken to- 
gether. Fig. 7 shows spinon creation upon a single spin ex- 
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Fig. 6. (color online) Representative temperature-dependent QMC results for bond ((a) and (d)), charge ((b) and (e)), and spin susceptibilities ((c) and (f)) 
for a 36 site zigzag ladder with U = 4, V = 0, and a = 0. /a/'l is 0.4 in panels (a)-(c), and 1.4 in panels (d)-(f). 



and high temperatures are dominated by high-spin configu- 
rations containing multiple spinons. Thus the signatures of 
spinons, including spinon binding can also are to be be found 
in the finite-temperature QMC results of Fig. 6. In the ID 
hmit spinons have local 4kf charge or bond order. As in 
the ID chain (reference^^), when f2/fi is small we find that 
the 4kf bond and charge susceptibilities increase with tem- 
perature (see Fig. 6(a)-(b) at <7 = n). 

When f2/fi is of order unity, several diff'erences as seen in 
the susceptibilities that can be correlated with the binding of 
spinons shown in Fig. 7. First, at intermediate temperatures, 
the charge and spin susceptibilities have a broad plateau at 
q a: 37r/4. When spinons bind, the larger size of these de- 
fects moves the charge and spin response away from 4kf and 
towards a smaller wavenumber Second, in this parameter re- 
gion the 4kf bond susceptibility remains small regardless of 
temperature and instead at higher temperatures a broad re- 
sponse from < q < n/2 appears (Fig. 6(d)). This small- 
q (long wavelength) bond order response is a result of in- 
creasing numbers of bound spinons created in high-spin con- 
figurations. Each high-spin configuration has multiple bound 
spinons equally separated from each other, giving a long- 
wavelength bond order distortion with Q < q < njl. Related 
to this is the observation ih^Axsi^kf) varies only weakly with 
temperature in the ladder limit. In the ID limit, each spinon 
separates regions where the PEC is out of phase by two lat- 
tice units (Fig. 7(a)); due to this phase difference the intro- 
duction of spinons will lead to a rapid decrease in the 2kf 
response. On the other hand in the ladder limit, the PEC re- 
mains "in phase" on either side of bound spinons (Fig. 7(b)), 
resulting in a weaker temperature dependence of the 2kf bond 
susceptibility — ^bound spinons suppress the BCDW locally, 
but do not disturb the overall phase of the density wave as 
do single spinons. 

5. Discussion 

We have shown that the | -filled coiTelated zigzag ladder 




Fig. 7. (color online) MPS-QMC charge and spin density in the triplet state 
versus site index j along the t\ direction. Parameters are W = 36, U = 6, 
V = I, and a = 1 .6. In panels (a) and (b), t2 / fi =0.2 and 1 .4, respectively. Cir- 
cles (diamonds) are chai'ge (spin) density. Arrows indicate location of spinon 
defects (see text). Statistical errors are smaller than symbol sizes. 



this lattice, the WC state driven by strong nearest-neighbor 
Coulomb interactions in ID is strongly suppressed. Instead, a 
single PEC state with CO pattern ■ ■ ■ 1 100 ■ ■ ■ occurs over a 
wide range of fa/fi (0 < fa/fi < 1.707). Second, unlike in the 
ID i-filled band the bond order pattern in the zigzag ladder 
is S MWM for all Coulomb interactions, and the bond distor- 
tion SWSW never occurs. Third, with increasing f2/fi from 
the ID limit, the CO amplitude decreases for fixed V until it 
vanishes above a critical value. Surprisingly, the singlet-triplet 
gap increases in magnitude at the same time. We have shown 
that this increase is accompanied by the binding of spinon ex- 
citations. The singlet-triplet gap can be very large when f2/fi 
and V are both moderately large (see Fig. 5(d)). Finally, our 
QMC calculations show that in the zigzag ladder there is no 
distinct intermediate state between the low temperature PEC 
and the high temperature metallic state, and there is a single 
metal-insulator transition that is accompanied by simultane- 
ous bond and charge distortions and spin gap. The existence 
of simultaneous transitions is expected here as the spin-singlet 
here, is of the "intra-Himer" tvne as shown in Fip_2(c)-(rlV 
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In the isotropic ^-filled band 2D square lattice, adding a 
frustrating diagonal /' bond results in a transition from a uni- 
form state to a paired PEC state once t! exceeds a critical 
value. As with spin systems, the zigzag ladder provides a 
simple model for the study of the frustration-driven PEC state 
in 2D, with the difference that the frustration does not create 
the SG state in the ladder, but rather enhances the SG that is 
already present in the unfrustrated model. This enhancement 
can be large due to the cooperative interaction between the 
two kinds of electron-phonon interactions that are possible in 
the ladder, as shown in Fig. 3. We discuss below the impli- 
cations of our work for understanding experiments on | -filled 
materials in general (including both 2D and 3D). We first con- 
sider specifically those materials which have been suggested 
to be ladders based on their crystal structures. We then also 
consider the broader class of ^-filled band materials. Ideally, 
this second class of materials require understanding of the ex- 
citations and thermodynamics of the higher dimensional PEC, 
which is a much more formidable task than the ladder calcu- 
lations. We nevertheless point out that broad conclusions can 
be drawn for these systems based on our current work. The 
majority of the materials we consider belong to the large fam- 
ily of low dimensional organic charge transfer solids (CTS), 
within which many examples of SG ground states are found 
with quasi-lD and 2D lattice.^' We also point out a few inor- 
ganic ^-filled materials that exhibit similar transitions. 

5.1 I -filled ladder candidates 

While most ladder materials found have been i-fiUed,^^ 
the family of CTS materials (DTTTF)2M(nmt)2, where M 
is a metal ion, are likely candidates for | -filled band lad- 
ders. In this series of compounds, the M=Au and M=Cu 
(for both the metal ion is diamagnetic) have been studied in 
most detail.-'''"''^ Structurally, these materials consist of pairs 
of DTTTF stacks, each with ^ -filled band of holes, sepa- 
rated by stacks of M(mnt)2 which are Mott-Hubbard semi- 
conductors with i-filled electron bands. It is then likely that 
the each pair of DTTTF stacks behaves as ladders. Note that 
the stack direction corresponds to the ?2 direction within our 
model, and hence these systems lie in the parameter regime 
ti>h. The very large^^ Tsg in (DTTTF)2Au(mnt)2 (-70 K) 
and (DTTTF)2Cu(mnt)2 (~90K), which are nearly an order of 
magnitude larger than the spin-Peierls transition temperatures 
in the ID |-filled systems,^** supports the conjecture that these 
systems are ladders. 

The MI and SG transitions are, however, distinct in these 
compounds, which would argue against the zigzag ladder pic- 
ture, at least its simplest version. For M=Au, a broad MI tran- 
sition occurs at Tmi ~ 220 K, followed by a decrease in the 
magnetic susceptibility at 70 K.^"* Below the MI transition, 
difl'use X-ray scattering at fo*/2 indicates dimerization along 
ti, but broad line widths suggest the dimerization order is 
not long-ranged.^"* The M=Cu salt is isostructural to the Au 
salt with slightly smaller lattice parameters due to the smaller 
metal ion. The MI transition for M=Cu occurs at 235 K, and 
unlike the Au salt is a sharp, second-order phase transition 
that is accompanied by doubling of the unit cell in the lad- 
der direction.-''' Changes in the optical properties at Tmi and 
Tsg for the two salts are also difi'erent.^^'^^ For M=Au, at Tmi 
symmetry breaking occurs along the rung direction (perpen- 



dicular to the DTTTF stacks), while at Tsg, symmetry break- 
ing is predominantly along the stacks. In contrast, optical 
response indicates symmetry breaking in both rung and stack 
directions at Tmi for M^Cu.-'^ 

We believe that the j -filled band zigzag ladder model is 
nevertheless a valid description for both M = Au and Cu 
at low temperatures. The only other competing model for 
these systems is the rectangular ladder model, wherein each 
DTTTF molecule is coupled to a single other such molecule 
on the neighboring stack. Such a description would be against 
the known crystal structures.^-' Furthermore, within the rect- 
angular ladder model, there needs to occur a high tempera- 
ture metal-insulator transition accompanied by in-phase bond 
dimerization, such that each dimer of DTTTF molecules has 
a single electron; the ladder after dimerization would be akin 
to rectangular spin ladder, which has SG for all interstack 
spin exchange. The two stacks need to be identical within 
the model and hence there is no symmetry breaking within 
the rectangular ladder scenario along the rung direction at 
any temperature. Neither is there any CO within the model, 
in contradiction to what is found in optical measurements.^^ 

There can be several different reasons why the high temper- 
ature MI transition occurs within the zigzag ladder scenario. 
First, muon spin rotation experiments suggest significant in- 
terladder coupling,^^ that has been ignored in our isolated lad- 
der model. Second, our model does not take into account the 
temperature-dependent lattice expansion that is common to 
CTS crystals. Note that lattice expansion will affect the inter- 
stack hopping f 1 much more strongly than the intrastack hop- 
ping f2- It is then conceivable that at high temperatures the 
interstack distance is large enough (and t\ is small enough) 
that ttjtx is greater than the critical value 1.707 and the sys- 
tems behave as independent chains. The lattice contracts at 
reduced temperatures, increasing t\ and reducing fa/fi, when 
the systems exhibit zigzag ladder behavior. This would re- 
quire hjti close to 1.7 in the experimental systems, which 
is indeed close to ratio of the calculated hopping integrals 
for the M = Au system.-'^ As seen in Fig. 4(d), large f2/fi 
would be in agreement with the unusually large SG seen in 
the (DTTTF)2M(mnt)2. There are however additional compli- 
cations. During the synthesis of the (DTTTF)2M(mnt)2 salt, 
the 1 : 1 salt (DTTTF)M(mnt)2 is also produced and crystals of 
the 1:2 salt must therefore be separated from this mixture for 
experiments.^^ Relative to other CTS, available experimen- 
tal data is thus more limited. Experimental determinations of 
the pattern of the CO below Tsg (and above, if any), and of 
the temperature-dependent lattice distortion are necessary for 
resolution of the above issues. 

5.2 General classification of \ -filled materials 

As discussed in Section 3.1, in ID the SG and MI tran- 
sitions are distinct when the ground state broken- symmetry 
state has the bond pattern SWSW', but are coupled together 
in a single transition when the bond pattern is SMWM as oc- 
curs in the zigzag ladder we have considered here. In ID, 
the strength of e-e interactions determines which bond pat- 
tern is favored. As SG transitions are found in a number of 
I -filled materials with 2D as well as 3D lattices, general- 
izations of these results to higher lattice dimensionality are 
of great interest. Our results in Section 4 for the zigzag lad- 
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Material 


D-n 


TsG 


TcB 


Ref. 






(K) 


(K) 




MEM(TCNQ)2 


lD-2 


17 


335 


42,43 


(TMTTF)2PF6 


lD-2 


19 


70 


38,44 


6'-(ET)2RbZn(SCN)4 


2D-2 


20 


195 


45-48 


EtMe3P[Pd(dmit)2]2 


2D-2 


25 


> 300 


49,50 


/3"-DODHT)2PF6 


2D-2 


40 


255 


51,52 


(DMe-DCNQI)2Ag 


lD-2 


80 


100 


53 


Q'-NaV205 


2D-1 


34 


34 


54 


a'-(ET)2l3 


2D-1 


135 


135 


55,56 


(BDTFP)2PF6(PhCl)o.5 


lD-1 


175 


175 


57-59 


CuIr2S4 


3D-1 


230 


230 


60,61 


(EDO-TTF)2PF6 


lD-1 


280 


280 


62,63 



Table I. | -filled materials with spin gapped ground states. For each the di- 
mensionality and number of transitions D-n, spin-gap transition temperature 
TsG, charge-bond ordering temperature Tcb, and references are listed. 



Strength does not necessarily determine the bond distortion 
pattern and thermodynamics — ^lattice structure and frustration 
are also important. 

We point out an empirical criterion here that rationalizes 
separate versus coupled SG-MI transitions in | -filled band 
materials at large, that we arrive at by simply extrapolating 
from the ID and zigzag ladder results. Our observation is that 
if the low temperature structure is such that the singlet bond 
is interdimer, and the strongest bond is between intradimer 
charge-rich and charge-poor sites (Fig. 2(b)), there occur dis- 
tinct transitions involving charge and spin degrees of freedom. 
Conversely, if the intradimer bond is between a pair of charge- 
rich sites (Fig. 2(c)), and the SG is due to intradimer spin- 
singlets, there is a single coupled SG-MI transition. In this 
latter case in general Tsc is high. The first of these two obser- 
vations was noted previously by Mori."*" We do not have any 
microscopic calculation to justify these conclusions for 2D 
and 3D; they are based on the mappings of Fig. 2. In Table I 
we give a hst of materials for which the bond and/or charge 
ordering pattern below Tsg is known. In all cases our sim- 
ple criterion appears to be valid (see below). Two of the en- 
tries require additional explanation. (TMTTF)2PF6 is already 
insulating at room temperature because of lattice dimeriza- 
tion; the high temperature CO at 70 K here is of the WC type, 
but there is strong redistribution of the charge below Tjc/' 
clearly placing this system into the category of materials that 
undergoes two transitions.^^ The low Tsg is also a signature 
of this. The situation with a'-NaV205 is exactly the opposite. 
This system is also insulating already at high temperature, but 
now the charge-ordering and spin-gap transitions occur at the 
same temperature, indicating that these transitions are cou- 
pled. While Table I is meant to be representative and not com- 
prehensive, we are unaware of examples where our criterion 
fails. We describe individual materials in detail below. 

5.2.1 Two transitions: inter-dimer singlet formation 

We have previously discussed quasi- ID CTS where two 
transitions occur.^^ MEM(TCNQ)2 is one example where the 
low temperature bond pattern has been well characterized by 
neutron scattering."'^ In MEM(TCNQ)2, the MI transition oc- 
curs near room temperature at Tmi=335 K, followed by a SG 
transition at Tsg=17.4 K.^^-"'^ The bond pattern at low temper- 



ature is SWSW'.^ The (TMTTF)2X series is another quasi- 

ID example where CO and SG occur at different tempera- 
tures. 44. 6«6 

In the 2D 0-(BEDT-TTF)2MM'(SCN)4 series, CO occurs 
at the high temperature MI transition'*^'46 followed by the SG 
transition at Tsg ~20 k.4^'4'^'48 X-ray studies in the tempera- 
ture range Tsg <T<Tco show that the strongest bond orders 
in the 2D BEDT-TTF layers are between molecules with large 
and small charge density (see Fig.7 in reference4^). These and 
lower temperature X-ray results'*** have indicated that the spin- 
singlets in the SG phase are located on the inter-dimer bonds, 
as in the SWSW' bond pattern in ID. 

Yet another 2D material that very clearly shows two tran- 
sitions and in which the charge-bond distortion patterns are 
known is EtMe3P[Pd(dmit)2]2.''^'^*' The material is semicon- 
ducting already at 300 K (Tmi > 300 K). The magnetic 
susceptibility at high temperatures corresponds to that of a 
Heisenberg antiferromagnet on a triangular lattice with J = 
250 K. Below a relatively low Tsg=25 K the system enters a 
distorted phase, with the intermolecular bond distortion pat- 
tern clearly of the dimerized dimer type, and the strongest 
bonds between charge-rich and charge-poor sites. This ma- 
terial undergoes superconducting transition under pressure.^^ 

5.2.2 Coupled transitions: intra-dimer singlet formation 

In cases where Tmi and Tsg coincide, the spin gap 
transition temperature tends to be quite high. Examples 
here include (EDO-TTF)2X, which shows a first order 
Ml transition at high temperature, 280 K for X=PF6 
and 268 K for X=AsF6.^2 This transition coincides with 
TsG-*^ Optical experiments determined that the charge or- 
der pattern in the low temperature phase is •••llOO---, 
with the strongest bond between molecules with large 
charge density.^^ A coupled SG-MI transition occurs at 
Tmi=175 K in (BDTFP)2PF6(PhCl)o.5." While structurally 
(BDTFP)2PF6(PhCl)o.5 appears to be ladder-Uke,'^'^'^ X-ray^^ 
and optical measurements^^ show that in the low-temperature 
phase tetramerization takes place along the stacks with the 
SMWM bond pattem.^^ The coupled SG-MI transitions 
found in these two materials are consistent with our criterion 
above. 

Beyond quasi- ID materials, in the 2D CTS examples can 
be found with coupled transitions, which also typically take 
place at a relatively high temperature. In a-(BEDT-TTF)2l3 
Tsg=135 K coinciding with the MI transition.^^ Similar to 
(EDO-TTF)2X, the SG-MI transition is first order and coin- 
cides with a large structural change.^^ In the low temperature 
phase, the strongest bond is again between the sites of largest 
charge density.^*" 

Similar transitions are observed in inorganic ^-filled ma- 
terials as well. The inorganic spinel CuIr2S4 is one example 
in which the Ir-ions form the active sites with | -filled elec- 
tron band (^-filled hole band).*'* In Culr2S4 a coupled SG-MI 
transition occurs at 230 K, below which the criss-cross chains 
of Ir-ions are charge-ordered as Ir4+-lr"*^-lr^^-Ir^^, with the 
strongest bonds between the spin i hole-rich Ir4+ ions and 
weakest bonds between the spin hole-poor Ir-'"^ ions..^"'^* A 
more complex coupled transition occurs in Q''-NaV205, where 
a coupled CO-SG transition occurs at 34 K within an insulat- 
ing phase. Structurally a'-NaV205 consists of rectangular V- 
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ion based ladders linked by zigzag V-V bonds. Below the tran- 
sition the V-ions are charge disproportionated and there oc- 
curs a period-4 V'^^-V'^^-V^^-V^^ CO within the zigzag links 
between the rectangular ladders. Once again, the strongest 
bonds are between the spin j electron-rich V"*^ ion pairs and 
the weakest bonds between the spin electron-poor V^^ ion 
pairs, in agreement with our criterion for coupled CO-SG 
transition. 

5.3 Possible relationship with superconductivity 

We have recently suggested that superconductivity in 
strongly-correlated | -filled systems is due to a transition from 
an insulating PEC state to a paired-electron liquid.''' '^ ™ 
Within this model, the spin-singlets of the PEC become mo- 
bile with further increase in frustration. The fundamental the- 
oretical picture is then analogous to bipolaron theories of su- 
perconductivity,^' with two diff'erences: (i) the pairing in our 
model is driven by antiferromagnetic correlations (as opposed 
to very strong e-p interactions that screen out the short-range 
Coulomb repulsion), and (ii) nearly all the carriers are in- 
volved in the pairing. The effective mass of the spin-bonded 
pairs is an important parameter, and overly strong binding 
will reduce pair mobility. This would suggest that supercon- 
ductivity is more likely in materials with inter-dimer singlets. 
As noted in the above, such systems tend to have the dimer- 
ized dimer structure. In contrast, in those materials with intra- 
dimer singlets that form at high temperature, the stronger pair 
binding would lead to a pair mobility too low to achieve 
superconductivity-in these cases the ground state would re- 
main an insulating spin-gapped PEC with charge and bond 
order. It is interesting to note that that such a correlation was 
suggested from empirical observations alone by Mori.'*" 
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Appendix A: MPS-QMC method 

Matrix-product states are best known as the class of wave- 
functions ultimately reached by the DMRG process^^ and 
provide highly eflicient representations of wavefunctions for 
quasi-lD quantum systems. An alternate set of methods to 
DMRG has recently emerged which combines a MPS wave- 
function representation with MC sampling.^''' In this ap- 
proach, rather than the renormalization procedure used by 
DMRG, Monte Carlo (MC) sampling is used to evaluate ex- 
pectation values of the MPS wavefunction. The use of MC 
has certain advantages, such as potentially better computa- 
tional scaling,^^ and the ease with which the method can be 
paralleUzed through trivial parallelization of the MC averag- 
ing. 

The MPS-QMC method we use is described in Refer- 
ence,^^ where it was applied to the ID quantum Ising model 
in a transverse field. Here we show that the method can suc- 
cessfully be applied to more complicated Fermion Hamiltoni- 




Fig. A-1. (color online) Relative error in the ground state energy using 
MPS-QMC as a function of number of QMC steps k and matrix size D for 
a 20 site ladder with periodic boundary conditions and t[ = t2 = i, U = 6, 
V = I, and a = 0. Calculations for D=4, 8, and 16 were single optimization 
runs, while those for 32 and 48 were restarted after 50-75 steps (see text). 



system is written as 

m = 2Tr[Ai(s,)A2(s2) ■ ■ ■ An(sn)]|siS2 ■ ■ ■ Sn). (A-1) 

In Eq. A-1, \S} - \s1s2 ■ ■ ■ sn) is a single many-electron con- 
figuration where the state on each site / is 5,. In the spin-^ 
model considered in Reference^^ i, = +1. For the Hubbard 
model Sj takes four possible values corresponding to either 
empty, spin-up, spin-down, or doubly occupied states. Aj are 
D X D matrices, where D is an adjustable parameter. As we 
study ground state configurations containing spinon defect 
states that break translational symmetry, we keep separate ma- 
trices Ai(si) for each lattice site. In the MPS-QMC method the 
matrix elements of the matrices A,(s,) are stochastically opti- 
mized to minimize the expectation value of H. 

The method may be summarized as follows: the matrix el- 
ements a'^jisit) are first initialized to hold random values. Ma- 
trices are normalized so that their Frobenious norm is unity, 
^JTr(AA'^) - 1. The expectation values of the energy and 
derivatives of the energy with respect to a'^j{sk), k - 1 ■ ■ - N, 
are evaluated using MC sampling and a Metropolis accep- 
tance probability. Configurations are sampled according to the 
weig ht W2(5) with 

W{S) = Tr[Ai(si)A2(s2) - - - An(sn)]- (A-2) 

The MC updates used to sample the configurations {S } inter- 
change electrons between neighboring sites, i.e. { 2) ^ {]^ 
). In such a process, two matrices in Eq. A-2 are changed. 
The changed weight, W(S'), may be efficiently generated by 
saving a series of "right" and "left" matrix products, and se- 
quentially attempting updates^*" for sites / = I ■ ■ - N. The es- 
timates of the derivatives are used to optimize the initial ran- 
dom matrices using a stochastic optimization scheme.^'' The 
calculation is divided into a series of steps = 0, 1, 2, - - - . At 
each step, the amplitude of the stochastic noise is decreased, 
and the number of MC samples is increased to provided in- 
creasingly accurate derivatives as the minimum energy is 
approached. The overall operation count of the MPS-QMC 
method scales more favorably with the matrix size (oc D^) for 
periodic systems than does traditional DMRG.^*" 

Fig. A-1 and Table A-1 compare the accuracy of MPS- 
QMC against exact diagonalization for the energy, kinetic en- 
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D 


E 


KE 


PE 




4 


-0.92942(1) 


-1.46896(4) 


0.53954(4) 


0.15115(3) 


8 


0.940809(7) 


-1.46732(4) 


0.52651(4) 


0.15280(3) 


16 


-0.942555(2) 


-1.46824(4) 


0.52568(4) 


0.15608(2) 


32 


-0.9426859(6) 


-1.46818(3) 


0.52549(3) 


0.15647(4) 


48 


-0.9427005(5) 


-1.46821(3) 


0.52551(3) 


0.15649(4) 


ED 


-0.9427135 


-1.468223 


0.525510 


0.156520 



Table A-1. Comparison of the ground state energy per site, kinetic energy 
per site, potential energy per site, and spin structure factor 5o-(?r/2) between 
MPS-QMC and exact diagonalization a 20 site ladder with periodic boundary 
conditions and f2/'i = 1, f = 6, and V = I. Statistical Monte Carlo sampling 
errors in the last digit are shown in parenthesis. 



ergy, potential energy, and spin structure factor, defined as 

Saiq) = ^ ^"^'"'^<(«iT - «ii)(«/T - «a))- (A-3) 



jl 



Fig. A l shows the relative error in the ground state energy 
as a function of the number of sampling blocks for a 20-site 
ladder with t\ - t2 - I, U - 6 and V = 1. Each step k of the 
algorithm is further divided into a number n of blocks consist- 
ing of m Monte Carlo updates each, with the matrices A,(5' ,) 
updated using the calculated derivatives after each block, n 
and m are increased as n = kn() and m - kmo, with no = 25 
and mo = 100 in Fig. A-1. We find that larger matrix sizes are 
needed here compared to those typically used for spin mod- 
els; this is expected because of the greater number of quan- 
tum states per site. However, even with the much larger num- 
ber of parameters (matrices are also site-dependent in our cal- 
culations) the stochastic optimization method performs well. 
Comparing with N - 20 site exact results, relative energy ac- 
curacy of order 10"^ is achievable with D = 48. As shown in 
Table A l, local quantities such as the kinetic energy (bond 
order) and potential energy converge quite rapidly with D. As 
expected, long-range correlations require larger D for simi- 
lar numerical accuracy. The 36-site calculations presented in 
Section 4 used D = 48 matrices. 

In order to treat the e-p terms in Eq. 1, we use the an it- 
erative self-consistent method^^ to update the bond distor- 
tions Aij. The self-consistency equations for the lattice are 
updated once every MC block following the calculation of 
charge densities and bond orders. Self-consistent MPS-QMC 
results were verified against exact self-consistent calculations. 

Appendix B: Finite-size scaling 

For the results reported in Figs. 4 and 5 we performed 
finite-size scaling from calculations on 20, 28, 36, and 60 site 
ladders. The 20 site data was from exact lanczos calculations 
and larger system data from the MPS-QMC method detailed 
in Appendix A using a matrix size of D = 48. Here we give 
further details of our extrapolation procedure. 

Because the properties of quasi-lD systems display an al- 
ternation with as the chain length increases, we have chosen 
only ladders with An + 2 electrons. These systems also have 
closed shells in the noninteracting limit. The PEC ground 
state further requires that be a multiple of four Fig. B-l(a) 
shows the finite-size extrapolation for the charge dispropor- 
tionation An for both small and large f2/fi. An as well as r4 
and Es (not shown here) change only by a small amount from 
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Fig. B-1. (color online) Finite-size scaling of observables for U = 6, V = 
1, and a = 1.6. Circles (squares) are for tilti = 0.2 (b/'i = 1-6). N is the 
number of lattice sites. Filled symbols show the extrapolated values of the 
observables. Panel (a) shows the charge disproportionation Ah and panel (b) 
the singlet-triplet gap Ast- Lines are least-squares fits to the points. 



singlet-triplet gap Ast- Errors in the extrapolation of Ast are 
larger than those in An, r4, and E^, because the gap is calcu- 
lated from the energy difference of separate 5=0 and 5 = 1 
calculations, each of which have different lattice configura- 
tions that are optimized self-consistently. This accounts for 
the larger error bars in Fig. 4(d) compared to Fig. 4(a)-(c). 
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